function  [GG, RR, CONS, eu, SDX, ZZ, initss, ssvec, flag, ssnames , stanames ,shonames]=...
    modchimodstama(estpr,options,op,indp,calpr,con)

% =========================================================================
% CHIMOD ST AMA 
% AMA version of that model 
% Single trend model, with a single price observable and real Y,C,I and W 
% using that same price series. 
% Formerly AYPISTARWNX 
% Policy responds to annualized output growth and annualized inflation. 
% There is a split in the labor disutility shocks, and the AR part feeds
% into X 
% 
% Decouple LSS and LSScale==1
% Alejandro Justiniano March 16 2010 
% 
% Checked solution is IDENTICAL to within 1e-12 to @MODCHIMODST
% =========================================================================
if nargout < 11 
    flag_repstanames=0; 
else 
    flag_repstanames=1; 
end 
if nargout < 10 
    flag_repssnames=0; 
else
    flag_repssnames=1; 
end

flag.ssok=1;
flag.solver=1; 

initss=[]; 
ssnames=[]; 

% INDEX for Ys
% ------------
y=1;            % output
k=2;            % capital
L=3;            % hours
Rk=4;           % return on capital
w=5;            % real wages
p=6;            % inflation
s=7;            % marginal cost
lambda=8;       % multiplier
c=9;            % consumption
R=10;           % interst rate
u=11;           % capital utilization
phi=12;         % multiplier
i=13;           % investment
kbar=14;        % "gross" capital
x=15;           % marginal utility of labor
z=16;           % productivity shock
g=17;           % public spending shock
miu=18;         % investment shock
lambdap=19;     % mark-up shock
psi=20;         % preference shock
b=21;           % discount factor shock 
% Lagged Variables 
z_1=22; 
z_2=23; 
y_1=24; 
y_2=25; 
y_3=26; 
p_1=27; 
p_2=28; 
% Observables
yobs=29;     
cobs=30;
iobs=31;
wobs=32;
psiwn=33;
mp=34;
dpma=35; 
pitarg=36;         % preference shock

NY=36; 
NETA=0; 

% ========================================
% INDEX for Shocks  
% ========================================
NX=9; 
Rs=1;           % MP
zs=2;           % Technology 
gs=3;           % Public spending
mius=4;         % Investment     
lambdaps=5;     % Mark-up
psis=6;         % Leisure preference AR
bs=7;           % Discount factor
psiwnsh=8;      % Leisure preference WN
pitargsh=9; 

%=======================
% cexstar=8;
% lambdaexstar=9;
% phiexstar=10;
% Rkexstar=11;
% iexstar=12;
%========================
% Index for the parameters 
% -----------------------
ncof=        29;  
numpar=ncof+NX ;  % Number of parameters 

if nargin < 4 || isempty( calpr ) 
    param=estpr  ; 
else
    param=zeros(numpar,1); 
    param(indp)=estpr; 
    param(con)=calpr; 
end; 
% Non SD parameters 
% -------------------
alpha=param(1);     % share of capital in the prod. function
delta=param(2);     % capital depreciation rate
iotap=param(3);     % price indexation (=0 is static indexation, =1 is dynamic)
iotaw=param(4);     % wages indexation
gamma100=param(5);  % steady state growth rate of technology
h=param(6);         % habit formation parameter
lambdapss=param(7); % steady state of mark-up shock (pins down steady state of wages)
lambdaw=param(8);   % wage mark-up
%psiss=param(9);    % steady state leisure preference shock
Lss=param(9);       % steady state for leisure (the ss for leisure is pinned down by psiss. convenient to parameterize the model in terms of Lss)
pss100=param(10);   % steady state quarterly inflation (multiplied by 100)
%%%rss100=param(11);   % steady state quarterly real rate of interest (multiplied by 100)
Fbeta=param(11);    % weird parameter of SW that is a function of beta and the steady state quarterly real rate of interest (ensures that beta is less than 1)
gss=param(12);      % steady state of the public spending shock
niu=param(13);      % curvature of the utility function for leisure
xip=param(14);      % price stickiness
xiw=param(15);      % wage stickiness
chi=param(16);      % elasticity of the capital utilization cost function
S=param(17);        % investment adjustment cost
fp=param(18);       % reaction to inflation
fy=param(19);       % reaction to output gap
rhoR=param(20); rhoz=param(21); rhog=param(22); rhomiu=param(23); rholambdap=param(24); rhopsi=param(25); rhob=param(26);
rhoARMAlambdap=param(27); 
if rhoARMAlambdap~=0; error('RHOARMALAMBDAP is not zero'); end 
rhopitarg=param(28); 
rhomp=param(29); %autocorrelation of MP shock


% ----------------------------------------------------------------
% Standard deviations (no need to specify them in the new version)
% ----------------------------------------------------------------
% sdR=param(27); sdz=param(28); sdg=param(29); sdmiu=param(30); sdlambdap=param(31); sdpsi=param(32); sdb=param(33);


% STEADY STATE
% ----------------------------
Lsscale=1; 

gamma=gamma100/100;
beta=100/(Fbeta+100);
rss=exp(gamma)/beta-1;
rss100=rss*100;
pss=pss100/100;
gss=1/(1-gss);

Rkss=exp(gamma)/beta-1+delta;
sss=1/(1+lambdapss);
wss=(sss*((1-alpha)^(1-alpha))/((alpha^(-alpha))*Rkss^alpha))^(1/(1-alpha));
kLss=(wss/Rkss)*alpha/(1-alpha);
FLss=(kLss^alpha-Rkss*kLss-wss);
yLss=kLss^alpha-FLss;
kss=kLss*Lsscale;
iss=(1-(1-delta)*exp(-gamma))*kss*exp(gamma);
F=FLss*Lsscale;
yss=yLss*Lsscale;
css=yss/gss-iss;

% dispaj('Consumption share: ',css/yss ); 
% dispaj('Investment  share: ',iss/yss );
% dispaj('Capital     share: ',kss/yss ); 
% dispaj('Inflation   SS   : ',pss100  ); 
% dispaj('Real Rate   SS   : ',rss100  ); 
% dispaj('Observable Rss   : ',rss100+pss100); 
% dispaj('Rk          SS   : ',100*Rkss ); 
% dispaj('beta        SS   : ',beta ); 

ssvec  =[css/yss;iss/yss;kss/yss;rss100;100*Rkss]; 
if  flag_repssnames==1;
    ssnames={'c/y'  ;'i/y'  ;'k/y';  'rss100';'Rkss100'};
end
% System Matrices 
GAM0 =zeros(NY ,NY) ; 
GAM1 = zeros(NY,NY) ;
GAM2 = zeros(NY,NY) ;
PSI = zeros(NY,NX) ;
C = zeros(NY,1) ;

% eq 1, prod. function (y)
% ----------------------------------------
GAM1(y,y)=1;
GAM1(y,k)=-((yss+F)/yss)*alpha;
GAM1(y,L)=-((yss+F)/yss)*(1-alpha);

% ===
% GAM1(ystar,ystar)=1;
% GAM1(ystar,kstar)=-((yss+F)/yss)*alpha;
% GAM1(ystar,Lstar)=-((yss+F)/yss)*(1-alpha);
%===

% eq 2, cost minimization (L)
% ------------------------------------------
GAM1(L,Rk)=1;
GAM1(L,w)=-1;
GAM1(L,k)=1;
GAM1(L,L)=-1;

% ===
% GAM1(Lstar,Rkstar)=1;
% GAM1(Lstar,wstar)=-1;
% GAM1(Lstar,kstar)=1;
% GAM1(Lstar,Lstar)=-1;
% ===

% eq 4, Phillips Curve (p)
% ------------------------------------------
GAM1(p,p)=1;
GAM0(p,p)=-beta/(1+iotap*beta);
GAM2(p,p)=iotap/(1+iotap*beta);
GAM1(p,s)=-((1-beta*xip)*(1-xip)/((1+iotap*beta)*xip));
GAM1(p,lambdap)=-1;        % this coefficient is normalized to rescale the sd to an interpretable number
                           % originally was = to -((1-beta*xip)*(1-xip)/((1+iotap*beta)*xip))*lambdapss/(1+lambdapss);
                           % notice that this is the only equation where lambdap enters
                           
%===
%GAM1(Rstar,sstar)=1;
%===

% eq 5, marginal cost (s)
% --------------------------------------------
GAM1(s,s)=1;
GAM1(s,Rk)=-alpha;
GAM1(s,w)=-(1-alpha);

%===
% GAM1(sstar,sstar)=1;
% GAM1(sstar,Rkstar)=-alpha;
% GAM1(sstar,wstar)=-(1-alpha);
%===

% eq 7, consumption foc (c)
% --------------------------------------------
expg=exp(gamma);
GAM1(c,lambda)=(expg-h*beta)*(expg-h);
GAM1(c,b)=-(expg-h*beta*rhob)*(expg-h)/[(1-rhob)*(expg-h*beta*rhob)*(expg-h)/(expg*h+expg^2+beta*h^2)]; 
GAM1(c,z)=-(beta*h*expg*rhoz-h*expg);
GAM1(c,c)=(expg^2+beta*h^2);
GAM0(c,c)=-beta*h*expg;
GAM2(c,c)=expg*h;

%===
% expg=exp(gamma);
% GAM1(cstar,lambdastar)=(expg-h*beta)*(expg-h);
% GAM1(cstar,b)=-(expg-h*beta*rhob)*(expg-h)/[(1-rhob)*(expg-h*beta*rhob)*(expg-h)/(expg*h+expg^2+beta*h^2)];  
% GAM1(cstar,z)=-(beta*h*expg*rhoz-h*expg);
% GAM1(cstar,cstar)=(expg^2+beta*h^2);
% GAM1(cstar,ecstar)=-beta*h*expg;
% GAM2(cstar,cstar)=expg*h;
%===

% eq 8, euler equation (lambda)
% --------------------------------------------
GAM1(lambda,lambda)=1;
GAM1(lambda,R)=-1;
GAM0(lambda,lambda)=-1;
GAM0(lambda,p)=1;
GAM1(lambda,z)=rhoz;

%===
% GAM1(lambdastar,lambdastar)=1;
% GAM1(lambdastar,Rstar)=-1;
% GAM1(lambdastar,elambdastar)=-1;
% GAM1(lambdastar,z)=rhoz;
%===

% eq 9, capital utilization foc (Rk)
% -------------------------------------------
GAM1(Rk,Rk)=1;
GAM1(Rk,u)=-chi;

%===
% GAM1(Rkstar,Rkstar)=1;
% GAM1(Rkstar,ustar)=-chi;
%===

% eq 10, capital foc (phi)
% ------------------------------------------
GAM1(phi,phi)=1;
GAM0(phi,phi)=-beta*exp(-gamma)*(1-delta);
GAM1(phi,z)=rhoz;
GAM0(phi,lambda)=-(1-beta*exp(-gamma)*(1-delta));
GAM0(phi,Rk)=-(1-beta*exp(-gamma)*(1-delta));

%===
% GAM1(phistar,phistar)=1;
% GAM1(phistar,ephistar)=-beta*exp(-gamma)*(1-delta);
% GAM1(phistar,z)=rhoz;
% GAM1(phistar,elambdastar)=-(1-beta*exp(-gamma)*(1-delta));
% GAM1(phistar,eRkstar)=-(1-beta*exp(-gamma)*(1-delta));
%===

% eq 11, investment foc (i)
% -------------------------------------------
GAM1(i,lambda)=1/(S*expg^2);
GAM1(i,phi)=-1/(S*expg^2);
GAM1(i,miu)=-1/(S*expg^2);
GAM1(i,i)=(1+beta);
GAM1(i,z)=(1-beta*rhoz);
GAM0(i,i)=-beta;
GAM2(i,i)=1;

%===
% GAM1(istar,lambdastar)=1/(S*expg^2);
% GAM1(istar,phistar)=-1/(S*expg^2);
% GAM1(istar,miu)=-1/(S*expg^2);   
% GAM1(istar,istar)=(1+beta);
% GAM1(istar,z)=(1-beta*rhoz);
% GAM1(istar,eistar)=-beta;
% GAM2(istar,istar)=1;
%===

% eq 12, capital utilization (k)
% -------------------------------------------
GAM1(k,k)=1;
GAM1(k,u)=-1;
GAM1(k,z)=1;
GAM2(k,kbar)=1;

%===
% GAM1(kstar,kstar)=1;
% GAM1(kstar,ustar)=-1;
% GAM1(kstar,z)=1;
% GAM2(kstar,kbarstar)=1;
%===

% eq 13, capital accumulation (kbar)
% ------------------------------------------
GAM1(kbar,kbar)=1;
GAM1(kbar,miu)=-(1-(1-delta)*exp(-gamma));          
GAM1(kbar,i)=-(1-(1-delta)*exp(-gamma));
GAM2(kbar,kbar)=(1-delta)*exp(-gamma);
GAM1(kbar,z)=(1-delta)*exp(-gamma);

%===
% GAM1(kbarstar,kbarstar)=1;
% GAM1(kbarstar,miu)=-(1-(1-delta)*exp(-gamma));
% GAM1(kbarstar,istar)=-(1-(1-delta)*exp(-gamma));
% GAM2(kbarstar,kbarstar)=(1-delta)*exp(-gamma);
% GAM1(kbarstar,z)=(1-delta)*exp(-gamma);
%===

% eq 14, goods market clearing (u)
% ------------------------------------------
GAM1(u,c)=css/yss;
GAM1(u,i)=iss/yss;
GAM1(u,y)=-1/gss;
GAM1(u,g)=1/gss;
GAM1(u,u)=kss*Rkss/yss;

%===
% GAM1(ustar,cstar)=css/yss;
% GAM1(ustar,istar)=iss/yss;
% GAM1(ustar,ystar)=-1/gss;
% GAM1(ustar,g)=1/gss;
% GAM1(ustar,ustar)=kss*Rkss/yss;
%===

% ========================================
% OLD way of writing the wage PC
% Normalized for Unit Coefficient in wages 
% ========================================
% eq 15, wage curve (wtilda)
% AUX=(1+niu*(1+1/lambdaw))/(1-xiw*beta);
% GAM1(wtilda,wtilda)=1;
% GAM1(wtilda,x)=-1/AUX;
% GAM1(wtilda,p)=iotaw*xiw*beta;
% GAM1(wtilda,z)=-xiw*beta*(rhoz-iotaw);
% GAM1(wtilda,ewtilda)=-xiw*beta;
% GAM1(wtilda,ep)=-xiw*beta;
% % eq 16, MU labor (x)
% GAM1(x,x)=1;
% GAM1(x,psi)=-AUX*(1+beta)*xiw/(1-xiw);      
% GAM1(x,b)=-1/[(1-rhob)*(expg-h*beta*rhob)*(expg-h)/(expg*h+expg^2+beta*h^2)];
% GAM1(x,L)=-niu;
% GAM1(x,lambda)=1;
% GAM1(x,w)=-niu*(1+1/lambdaw);
% % eq 17, real wages indexation (w)
% GAM1(w,w)=1;
% GAM1(w,wtilda)=-(1-xiw);
% GAM1(w,p)=xiw;
% GAM1(w,z)=xiw;
% GAM2(w,w)=xiw;
% GAM2(w,p)=iotaw*xiw;
% GAM2(w,z)=iotaw*xiw;

% =========================================================================
% New way, also with unit coefficient 
GAM1(x,x)=1;
GAM1(x,w)=-1;
GAM1(x,lambda)=-1;
GAM1(x,b)=1/[(1-rhob)*(expg-h*beta*rhob)*(expg-h)/(expg*h+expg^2+beta*h^2)];
GAM1(x,L)=niu;
GAM1(x,psi)=1; 

% eq 16, wage (w): wage PC Written as in SW 
% -------------------------------------------------------------------------
slopewages=(1-beta*xiw)*(1-xiw)/((1+beta)*xiw*(1+niu*(1+1/lambdaw)));
GAM1(w,w)=1;
GAM0(w,w)= -beta/(1+beta);
GAM0(w,p)= -beta/(1+beta);
GAM1(w,x)=slopewages;
GAM1(w,z)=(1+iotaw*beta-rhoz*beta)/(1+beta);
GAM1(w,p)=(1+iotaw*beta)/(1+beta);
GAM2(w,p)=iotaw/(1+beta);
GAM2(w,w)=1/(1+beta);
GAM2(w,z)=iotaw/(1+beta);

% Removed 6/3/08 so models run before it should restore this normalization  
%GAM1(w,psi)  =   -1/4;  % Arbitrary normalization to faciliate a prior May 20 2008 
GAM1(w,psiwn)=     -1;

% =========================================================================

% eq 18, MP rule (R)
% ----------------------------------------
GAM1(R,R)=1;
GAM2(R,R)=rhoR;
GAM1(R,mp)= -1; 

GAM1(R,dpma)  = -(1-rhoR)*fp;
GAM1(R,pitarg)=  (1-rhoR)*fp;

GAM1(R,y)=-(1-rhoR)*fy/4;
GAM2(R,y_3)=-(1-rhoR)*fy/4;
GAM1(R,z)=-(1-rhoR)*fy/4;
GAM1(R,z_1)=-(1-rhoR)*fy/4;
GAM1(R,z_2)=-(1-rhoR)*fy/4;
GAM2(R,z_2)=(1-rhoR)*fy/4;

% GAM1(R,ystar)=(1-rhoR)*fy+fdy;
% GAM2(R,y)=-fdy;
% GAM2(R,ystar)=fdy;

% eq 19 - 24, exogenous shocks 
% ------------------------------------------
GAM1(z,z)=1; GAM2(z,z)=rhoz; PSI(z,zs)=1;
GAM1(g,g)=1; GAM2(g,g)=rhog; PSI(g,gs)=1;
GAM1(psi,psi)=1; GAM2(psi,psi)=rhopsi; PSI(psi,psis)=1;
GAM1(psiwn,psiwn)=1; PSI(psiwn,psiwnsh)=1; 
%%GAM1(ARMApsi,ARMApsi)=1; 
GAM1(miu,miu)=1; GAM2(miu,miu)=rhomiu; PSI(miu,mius)=1;
GAM1(lambdap,lambdap)=1; GAM2(lambdap,lambdap)=rholambdap; PSI(lambdap,lambdaps)=1;
GAM1(b,b)=1; GAM2(b,b)=rhob; PSI(b,bs)=1;
GAM1(mp,mp)=1; GAM2(mp,mp)=rhomp; PSI(mp,Rs)=1;
GAM1(pitarg,pitarg)=1; GAM2(pitarg,pitarg)=rhopitarg; PSI(pitarg,pitargsh)=1;

% eq 25 - 32, expectational terms
% ------------------------------------------

%===
% GAM1(ecstar,cstar)=1; GAM2(ecstar,ecstar)=1; PPI(ecstar,cexstar)=1;
% GAM1(elambdastar,lambdastar)=1; GAM2(elambdastar,elambdastar)=1; PPI(elambdastar,lambdaexstar)=1;
% GAM1(ephistar,phistar)=1; GAM2(ephistar,ephistar)=1; PPI(ephistar,phiexstar)=1;
% GAM1(eRkstar,Rkstar)=1; GAM2(eRkstar,eRkstar)=1; PPI(eRkstar,Rkexstar)=1;
% GAM1(eistar,istar)=1; GAM2(eistar,eistar)=1; PPI(eistar,iexstar)=1;
%===
GAM1(dpma,dpma)= -1;
GAM1(dpma,p)=1/4;
GAM1(dpma,p_1)=1/4;
GAM1(dpma,p_2)=1/4;
GAM2(dpma,p_2)= -1/4;

GAM1(y_1,y_1)=1; GAM2(y_1,y)=1;
GAM1(p_1,p_1)=1; GAM2(p_1,p)=1;
GAM1(p_2,p_2)=1; GAM2(p_2,p_1)=1;
GAM1(y_2,y_2)=1; GAM2(y_2,y_1)=1;
GAM1(y_3,y_3)=1; GAM2(y_3,y_2)=1;
GAM1(z_1,z_1)=1; GAM2(z_1,z)=1;
GAM1(z_2,z_2)=1; GAM2(z_2,z_1)=1;

% eq 33 - 36, lagged variables for the state-space repres.
% -----------------------------------------
% Observable output 

% yo(t)=gdp(t)-gdp(t-1)+z(t) 
% yo(t)-gdp(t)-z(t)= -gdp(t-1)
% gdp(t)=y(t)-( kss*Rkss/yss )*u(t) 
% yo(t) - y(t) + ( kss*Rkss/yss )*u(t)  - z(t) = -y(t-1) + (kss*Rkss/yss)*u(t-1) 
GAM1(yobs,yobs)=1;


GAM1(yobs,y)= -1; 
GAM1(yobs,z)= -1; 
GAM1(yobs,u)=kss*Rkss/yss; 

GAM2(yobs,y)= -1; 
GAM2(yobs,u)=kss*Rkss/yss; 

% Consumption 
GAM1(cobs,cobs)=1;
GAM1(cobs,c)= -1; 
GAM1(cobs,z)= -1; 

GAM2(cobs,c)= -1; 

% Investment 
GAM1(iobs,iobs)=1;
GAM1(iobs,i)= -1; 
GAM1(iobs,z)= -1; 

GAM2(iobs,i)= -1; 

% Real Wages 
GAM1(wobs,wobs)=1;
GAM1(wobs,w)= -1; 
GAM1(wobs,z)= -1; 

GAM2(wobs,w)= -1; 

% Standard Deviation 
SDX(1:NX,1:NX)=diag(param(ncof+1:end));
[GG, cc, RR, fmat, fwt, ywt, gev, eu] = amasolve(GAM0, GAM1, GAM2, zeros(size(GAM0, 1), 1), PSI);
flag.euok=isequal(eu,[1;1]);
% ================================================
% Matrix maps endogenous variables to observables 
% ================================================
ZZ=zeros(7,size(GAM0,2)); 
ZZ(1,yobs)=1; 
ZZ(2,cobs)=1; 
ZZ(3,iobs)=1; 
ZZ(4,L)=1; 
ZZ(5,wobs)=1; 
ZZ(6,p)=1;
ZZ(7,R)=1;
% ==========================
% Addition of the constants 
%===========================
% Note that they will be multiplied by zmat hence make sure to match the
% scaling 
CONS=zeros(NY,1); 
CONS(yobs)=gamma100;
CONS(cobs)=gamma100;
CONS(iobs)=gamma100;
CONS(L)=Lss;
CONS(wobs)=gamma100;
CONS(p)=pss100;
CONS(R)=pss100+rss100;

if  flag_repstanames==1;    
    stanames={'yhatwu','k','hours','Rk','what','dp','mgc','lambda','chat ','nomr','capu',...
        'phi','ihat','kbar','wdmrs','neutral','g','miu','pmark','psi','b',...
        'z_1','z_2','y_1','y_2','y_3','p_1','p_2',...
        'output','consumption','investment','rwages','wmark','mp','dpma','pitarg'};    
    shonames={'mp','neutral','g','miu','pmark','psi','b','wmark','pitarg'};    
end       